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y—i Abstract 

We present a symbolic-numeric approach for the anal- 
ysis of a given set of noisy data, represented as a finite 
<^ set X of limited precision points. Starting from X and a 

permitted tolerance e on its coordinates, our method auto- 
matically determines a low degree monic polynomial whose 
associated variety passes close to each point of X by less 
than the given tolerance e. 
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1 Introduction 

oo 
O 

It is a well-known matter that the analysis and modeling of real-world phe- 
nomena often relies upon measurements and observations which give rise to 
sets of data represented as affine sets of points. An algebraic way that exploits 
^ the collected data to construct a mathematical model of the observed phe- 

nomenon consists in computing the vanishing ideal of the affine set of points, 
that is the ideal comprising all polynomials which vanish at the given points. 
The vanishing ideal is classically determined using the Buchberger-Moller 
(BM) Algorithm [5J, a low complexity method which returns a Grobner basis 
of it. Nevertheless, the fact that data and relations from the real world are 
essentially characterized by a limited accuracy makes this algebraic approach 
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unfeasible in practice. In fact, even small perturbations of the points end up 
with very different vanishing ideals, which may have completely different 
bases, since the question whether or not a polynomial vanishes at a given set 
of points is obviously very sensitive to perturbations. For instance, this is a 
well-known phenomenon when using Grobner basis theory ([17). [IE]) which 
turns out to be unsuitable as a numerical tool. Thus, when dealing with real- 
world phenomena, the sole use of the classical algebraic techniques turns out 
to be insufficient, and a combination with numerical tools is required. 

Motivated by the need of describing real- world processes, the development 
of a numerical branch within a discrete mathematical discipline already took 
place in different areas of mathematics, such as Linear Algebra, Differential 
Equations, Optimization, etc. In the nineties, classical nonlinear algebra 
began to follow this trend, also thanks to the work of H. J. Stetter, whose 
book [21] is considered as a stepping stone at the interface of symbolic and 
numerical computation. This new emerging discipline has been given differ- 
ent names (see also [10]); in [22] its mathematical content has been clearly 
described as including the "areas of commutative algebra not over arbitrary 
fields or rings but over the analytically structured fields of the real or complex 
numbers", in which the real or complex metric is taken into consideration. 
Obviously, in the attempt of generalizing the notions of classical Commuta- 
tive Algebra to the real or complex coefficients, the greatest attention must 
be paid regarding the new concepts and terminology. In particular, in our 
case, the classical notion of vanishing ideal turns out to be inadequate for 
limited precision points, and so it needs to be replaced by an appropriate 
concept also suitable for numerical computations. 

In the literature a notion of vanishing ideal of a set of limited preci- 
sion points already exists; in [2], [10], [13], [15], [20], [23] it is defined in 
different ways as an ideal containing polynomials whose evaluations at the 
original points assume small values. Note that, according to this definition, 
the vanishing ideal of a set of limited precision points is not necessarily zero- 
dimensional; further, maximum relevance is given to the evaluations of the 
polynomials at the points, while the geometrical distance of the zero-locus of 
the polynomials from the original points is not taken into consideration. In 
this paper we follow a different approach: we give more relevance to the dis- 
tance points-variety and introduce the idea that the vanishing ideal of a set 
of limited precision points is generated by polynomials of low degree whose 
affine variety lies close to the original points. Among such polynomials, a 
crucial role is played by those of minimal degree, which provide an immedi- 



2 



ate interpretation of the phenomenon encoded by X and a simple geometrical 
representation of the data points, as shown in the following example. 

Example 1.1. Let X C M 2 be a set of points created by perturbing by less 
than 0.1 the coordinates of 10 points lying on the affine variety g = where 
g = - x - 2y + 2: 

X = {(0.95,1), (1.3, 0.5), (2.05, 1.98), (2.08,0), (3.18, -0.48), 

(5.05, 2.95), (5.05, -0.95), (7.2, -1.45), (9.98, 4), (10.05, -2)} 

Using standard techniques of Computational Algebra we obtain that the 
minimal degree of the polynomials vanishing at X is 4; this degree is too 
high to point out that the set X lies close to a parabola, while g and its 
zero-locus 2(g), almost crossing each point of X, provide a good numerical 
representation of the given data set. 

This paper, which is the first part of a wider investigation, reports on the 
problem of computing a polynomial / of low degree whose associated affine 
variety 2(f), which is defined as the set of points at which the polynomial / 
vanishes, almost crosses each element of the data set X. The final aim of 
our future work is to detect a set of polynomials whose minimal degree is 
strictly bounded by the minimal degree of the elements of the vanishing ideal 
of X and such that their zero-set is made up of points each differing from the 
corresponding point of X by less than the given tolerance. Put differently, 
we think to the numerical vanishing ideal as a truly zero-dimensional ideal 
whose exact zeros are small perturbations of the original limited precision 
points, and so numerically indistinguishable from them. An evident merit of 
this point of view is that the gain in simplicity due to the lower degree of the 
polynomials might offset the noise in the data and help to discover simpler 
laws that rule the phenomenon we aim to describe. 

In particular, in this paper we address the following problem: given a 
finite set X of points and a permitted positive tolerance e on the coordinates 
of each point, we look for a polynomial / whose degree is strictly bounded 
by the minimal degree of the elements of the vanishing ideal X(X) of X 
and whose affine variety 2(f) lies close to the points of X by less than e. 
We present a new algorithm (the Low-degree Polynomial Algorithm - LPA) 
that, starting from X and e, uses some peculiarities of the SOI algorithm j2] 
and NBM algorithm [10], both theoretically based on the BM algorithm, to 
determine a polynomial / of suitable degree. In most cases / also satisfies 
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some simple conditions, derived from Kantorovich theorem [23], under which 
the variety Z(f) is proved to lie close to the points of X by less than e. In 
these favourable cases / is a solution of our problem. 

More precisely, partly parallelizing the procedures presented in [2] and [TO] , 
we reformulate the addressed problem as the problem of solving an ordered 
finite sequence {Fi = 0}j of nonlinear systems subject to constraints. The 
first nonlinear system that admits an approximated solution satisfying the 
constraints allows the direct computation of the polynomial /. In the LPA 
we solve each nonlinear system Fi = using an iterative method which is a 
modified version of the classical Normal Flow algorithm [24]. Each iteration 
requires the solution of a linear system whose coefficient matrix is the eval- 
uation of the Jacobian Jp i at a suitable point. Since an analytic expression 
of Fi as well as of Jp i is computationally very hard to obtain, the evaluation 
of Jp i is performed using a new method that does not require the explicit 
knowledge of F^ Due to the technicalities of this method, its details are 
separately illustrated inside the Appendix. Moreover, since the evaluation of 
the Jacobian Jp i often turns out to be an ill-conditioned matrix, the corre- 
sponding linear system is solved using a method based on the Rank Revealing 
decomposition [14] . 

This paper is organized as follows. In order to be self-consistent and to fix 
notations, in Section [2] we introduce some basic concepts of Computational 
Commutative Algebra and Numerical Analysis. In Section[3]we introduce the 
modified version of the Normal Flow Algorithm to find an approximated so- 
lution of an underdetermined nonlinear system subject to constraints, paying 
particular attention to the ill-conditioned case. Section [4] presents the main 
results of the paper: the Low-degree Polynomial Algorithm, which performs 
the main steps of the method described above, and Theorem |4.3 which 
proves that, under certain hypotheses, the output of the LPA is a solution of 
our theoretical problem. In Section [5] we give numerical examples illustrating 
the functioning of the LPA. Finally, [6] contains the basic result for the evalu- 
ation of the Jacobian matrix associated with our nonlinear function, without 
passing through the explicit expression of it. 



2 Preliminary definitions and results 

We recall some basic notation about matrices. Let m, n > 1 be integers and 
Mat mxn (M) be the set ofmxn matrices with entries in K; if m = n we 
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simply write Mat n (R). Let A £ Mat mxn (M); we denote by ||A|| P the p-norm 
of A, p = 1, 2, oo. 

We divide the rest of this section into three subsections: the first one con- 
tains some preliminary definitions about the polynomial ring WL[x\, . . . ,x n ] 
and the notion of empirical points; the second one reports some basic re- 
sults on the numerical rank of a matrix; the third one describes a classical 
numerical approach to solve underdetermined nonlinear systems. 

2.1 Preliminary definitions 

We recall some concepts related to the polynomial ring P = M[xi, . . . ,x n ] 

(P!, HZ!). 



Definition 2.1. Let X = {pi, . 
of IR n , / be a polynomial and G 
polynomials. 



. . ,p s } be a non-empty finite set of points 
= {gi, . . . , g^} be a non-empty finite set of 



1. X(X) = {/ eP : f(pi) = Wpi £ X} is the vanishing ideal of X. 

2. Z(f) = {p £ M. n : f(p) = 0} is the affine variety associated with /. 

3. The IR-linear map evalx : P — > M s defined by evalx(/) = (f(pi), • • • , f(Ps)) 
is called the evaluation map associated with X. For brevity, we write 
/(X) to mean evalx (/) ■ 

4. The evaluation matrix of G associated with X, denoted by Mg(X) e 
Mat SX fc(M), is defined as having entry equal to gjipi), i-e. whose 
columns are the images of the polynomials gj under the evaluation map. 

Let t = xf 1 . . . x^ n , (3i e N, be a power product and O be a set of power 
products; we denote by d^t = PhX^ 1 . . . x^ k ~ l . . . x^ 1 the fc-th formal partial 
derivative of t and by dkO = {dkt : t £ O}. 

We formalize the idea of perturbed point by introducing a simplified ver- 
sion of the notions of empirical point, admissible perturbation and almost 
vanishing polynomials (see p], [2], [10], [21]). To this aim we recall that 
given a real value r\ 1 and a real function w(x) we say that w(x) = 0(rj k ), 
k £ N, if \w(x)\/r] k is bounded near the origin. 
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Definition 2.2. Let p,g6R™ be points and e > 0. 

1. The pair (p, e) = p e is called an empirical point: p is the specified 
value and e is the tolerance of p £ . 

2. The ^-neighborhood of p £ is defined as N(p £ ) = {p6R": ||p — < 
and contains all the admissible perturbations of p £ . 

3. Two empirical points p £ , q £ are said to be distinct if N(p £ ) nN(q £ ) = 0. 

Definition 2.3. Let X e = {pi, ■ ■ ■ ,p £ s } be a set of empirical points with 
uniform tolerance e and with X = {p\ , . . . ,p s } C M. n 

1. A set of points X = {p[, . . . ,p s } C W 1 is called an admissible per- 
turbation of X e if each p^ G N(p £ ). 

2. A polynomial g is called almost vanishing at X £ if ^(X)!^/!!*?!! = 0(e) 
where \\g\\ is the 2-norm of its coefficient vector. 

2.2 The numerical rank 

In this section we recall the definition of the numerical rank of a matrix 
A G Mat sxi (IR), and some basic results related to this topic ([H], [II])- In 
the case s > t such results will allow us to detect a partitioning of the columns 
of A (if s < t the partitioning involves the rows of A) into two submatrices 
Ai and A 2 (A 2 eventually the empty matrix) such that A x is well-conditioned 
and A 2 consists of columns (or rows) which are "almost" depending on the 
columns (or rows) of A x . 

Let A G Mat sx t(R); we denote by rank(A) the rank of A, by <Ji(A) the 
i-th singular value of A, by A< the pseudoinverse of A and by K 2 (A) = 
II ^4 II 2 1| ^4 1 1| 2 = cr i(^4)/ cr min{s,t}(^4) the 2-norm condition number of A. We 
denote by the vector or the matrix whose elements are equal to and 
whose dimension is deducible from the context. 

Given k > 1 and a small threshold 5 we say that A has numerical (S, k)- 
rank r if there exists a well-determined gap between a r (A) and a r+ \(A) and 
if both 5 and k5 lie in this gap. This concept can be formalized as follows. 

Definition 2.4. Let r < mm{s,t} and A G Mat sxt (lR); let 5 > and k > 1. 
The numerical (5, fc)-rank of A, denoted by rank^A), is equal to r iff 

o x > . . . > a r (A) > k5>5> a r+1 (A) >...> (r min{Stt} {A) 
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The singular values of A might not clearly split into small and large 
subsets making the determination of the numerical rank somewhat arbitrary. 
This leads to more complicated methods for estimating rank^A) which 
involve the LS problemJH]. In the following of the paper we assume that 
either the matrix A has full numerical rank or there exists a well-determined 
gap between cr r (A) and cr r+ i(^4). 

Remark 2.5. If A G Mat sx t(R) is a full numerical (5, fc)-rank matrix then it 
is well-conditioned, that is its condition number cannot be too elevated, as 
K 2 (A) = — ai ^ A \ A \ < ~Tjr^~ • On the contrary, if ranks J A) = r < minis, t] 
then the matrix is ill-conditioned, that is its condition number is elevated, 
since K 2 {A) = Jl(A > > 

The numerical (5, /c)-rank of A points out that there are no matrices with 
exact rank less than rank^A) which differ from A by less than S, as the 
following proposition shows. 

Proposition 2.6. Let A G Mat sxt (IR) with rank^A) = r. 

Then the set As = {M G Mat sx4 (lR) : \\M — A\\ 2 < 6} contains at least one 

matrix of rank r but no matrix of rank strictly less than r. 

Proof. From rank^A) = r it follows that a r (A) > 5 > a r+ i{A); further, 
from the Eckart- Young theorem [11] for each j = 1 . . . min{s, t} we have 
<7j(A) = min{||A - B\\ 2 : rank(5) = j - 1}. Let C G Mat sxt (R) be 
s.t. rank(C) = r and \\A — C\\ 2 = min{||A — B\\ 2 : rank(S) = r}; since 
a r+ i(A) < 5 then C G As. Further, since for each j < r we have Cj(A) > 5; 
it follows that As contains no matrix of rank strictly less than r. □ 

The numerical rank of A is preserved under small perturbations, since 
the following result about the sensitivity of the singular values holds. 

Proposition 2.7. Let A, E G Mat sx t(R); then for each 1 < j < min{s,t} 
we have \(jj{A + E) - <Tj(A)\ < <n(E) = \\E\\ 2 

Proof. See [II], Corollary 8.3.2 □ 



We end this section with Theorem 2.8 [14J which has important conse 



quences when applied to the case r = rank^A). 
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Theorem 2.8. Let r < t < s and A G Mat sX ^(R). There exists a permuta- 
tion matrix U G Matt(M) such that 

ATl=(Q 1 \Qi)( R > 1 j£) PI) 

where Q = (Qi \ Q 2 ) is orthonormal, Qi G Mat sxr (IR) ; Q2 £ Mat sx t_ r (IR) ; 
Rn G Mat r (IR) and R 22 £ Mati_ r (M) are upper triangular, and 

o-uun(Rn) > a j^l and \\R 22 \\ 2 < q(t, r)a r+l {A) (2.2) 
q(t,r) 

where q(t } r) = r(t — r) + min(r, t — r). 

Proof. See [H], Theorem 2.2. □ 



Theorem 2.8, when applied to the case r = rank^A), states the exis- 
tence of a set of r strongly independent columns of A or, equivalently, of 
a well-conditioned submatrix A\ G Mat sxr (IR) of A. In particular, it proves 
that there exist a permutation matrix II and a partitioning of AH = (Ai | A 2 ) 
where A 1 = QiRn G Mat sxr (R) and A 2 = QiR i2 + Q 2 R 22 e Mat sx t_ r (M), 
which has the following interesting numerical properties. If the gap between 
a r (A) and a r+ i(A) is large enough, that is if a r (A) ^> a r+ i(A), then A\ is a 
submatrix of A consisting of the maximum number of strongly independent 
columns w.r.t. the threshold S. In fact, from a r (Ai) = cr r (i?n) > > 0, 

we have that A± has full rank, and from Qi = A^R^ , H-R22II2 < q(t, r )°~r+i(A) 
and A 2 = Q\R\ 2 + Q 2 R 22 we have that each column of the matrix A 2 
can be expressed as the sum of a linear combination of the columns of 
A\ and a vector whose 2-norm is less than q(t, r)a r+ i(A). Furthermore, 



though A is ill-conditioned (see Remark | 2.5[ ), Ai is well-conditioned, since 



2.3 Underdetermined nonlinear systems 

Let D C R n , F : D — > IR' m , with m < n, be a different iable nonlinear 
function, Jf{x) be the Jacobian matrix of F at 2 and let F — be the 
underdetermined nonlinear system to solve. In this subsection we recall the 
Normal Flow Algorithm, a classical iterative method for approximating a 
solution of F = (see [21] and the references given there). 



8 



Algorithm 2.9. (The Normal Flow Algorithm - NFA) 

Let D C W 1 , F : D — >■ IR m be a differentiable nonlinear function, u < 1 
be a fixed threshold, and x E D be the initial point. Consider the following 
sequence of steps. 

NF1 Let h = (1,...,1)*. 

NF2 While \\h\\ 2 > uj 

• Compute the minimal 2-norm solution of Jp(x)h = —F(x). 

• Let x = x + h. 

Return x and stop. 

At each iteration the minimal 2-norm solution of the underdetermined 
linear system Jp{x)h = —F(x) can be either computed using Jp(x) or, more 
efficiently, performing a QR decomposition of Jf{x) (see [9]). 

The next classical result (see [21]) is a local convergence theorem for 
the NFA that also provides sufficient conditions for the existence of a local 
solution of F = 0. 

Theorem 2.10. (Kantorovich theorem) 

Let D C ]R n be an open set, F : D — > IR m be a differentiable nonlinear 
function and Jf{%) be of full rank m in an open convex set Q C D. Let 
xo G Q, 7] > and Q n = {x E Q : \\y — x\\2 < r] =3- y G Q}; suppose that 

• 3 7 > ; p G (0, 1] such that \\Jp(y) — J^(a;)||2 < -y\\y — x\^, Wx, y G Q 

• 3 fi > such that \\Jp(x)\\2 < fi, Wx E Q 

Then there exists (3 > satisfying r = 7M 1 ^ /3P < 1 and f~ < V su °h that if 
Xq E Q v and ||F(x )||2 < P then the iterates {xk}k=o,i,... determined by the 
NFA are well defined and converge to a point x* E Q such that F(x*) = 0. 

Proof. See HQ, Theorem 2.1. □ 

3 Underdetermined and ill-conditioned non- 
linear systems 

As reported in the Introduction, we reformulate the problem addressed in this 
paper as the problem of solving an ordered finite sequence of underdetermined 
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nonlinear systems subject to constraints which, in general, turn out to be ill- 
conditioned; in this section we present a method for solving such kind of 
nonlinear systems. 

Let D C W 1 be a set containing the origin, let F : D — > W n be a 
nonlinear function of class C 2 (D), and let Q £ = {x e R n : ||a?||oo < s}; 
our aim is to compute a solution of the underdetermined nonlinear system 
F = subject to x G Q £ fl D. As reported in Section 2J3 a classical method 
for finding an approximated solution of F = is given by the NFA, which 
returns successively better approximations of a zero of the system F = by 
computing at each step the minimal 2-norm solution of Jp{x)h = —F(x), 
where x is a suitable point of W 1 . Obviously, when the nonlinear system 
F = is ill-conditioned, that is when the condition number of the Jacobian 
matrix Jp of F evaluated in the constrained region is too elevated, the NFA 
could end up with unreliable solutions. In order to overcome the numerical 
instabilities in the computations and avoid the unreliable solutions which 
can occur in the ill-conditioned case, we present an alternative method that 
firstly replaces F with a new suitable function F slightly differing from F in 
the constrained region, that is such that — F(x)||2 is small on Q £ fl D, 

and well-conditioned at the origin (and also in a neighbourough of it). If 
this latter condition is not met, that is if the Jacobian matrix Jp(0) is ill- 
conditioned, we prove under a simple additional hypothesis that the original 
system F = subject to constraints has no solution (see Theorem 3.3). 
Successively our method applies a a modified version of the NFA to the 
system F = 0; the iterations of this modified version are stopped either when 
the computed deplacement is small enough, as in the classical version, or 
when the current iteration goes out of the region where Jp is well-conditioned. 
At the end this new algorithm returns either the current iterate x if the 
given constraints are satisfied, or a warning message. This method, that 
we present as the Root Finding Algorithm (see Algorithm 3.2), turns out 
to be backward stable [8] by construction, since it tries to solve, instead of 
the original constrained system, the nearby problem F = with the same 
constraints. 

In order to define the new function F : D — )■ IR m , we use the numerical 
rank of the matrix (Jf(0) | F(0)) as follows. Given a threshold 5 > s of the 
same order of magnitude as e and an integer k > 1, we detect the numerical 
(S, /c)-rank of the matrix (Jp(0) | F(0))* by means of the SVD decomposition. 



Let IT be the permutation matrix whose existence is stated in Theorem 2.8 
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Rn R12 
R 22 



(3.1) 



with r = rank,5 ]fc ( J F (0) | F(0))*, which satisfies 

f§£ ) n = Wl 1 g2 

where Q = (Qi \ Q 2 ) is orthonormal, Q 1 G Mat n+ i xr (IR), Q 2 G Mat„ + i xm _ r (IR 
Rn G Mat r (R) and R22 G Mat m _ r (IR) are upper triangular. Further we 
have cr miri (i?n) > (k/q(m,r))5 and H-R22II2 < q(m,r)5, where q(m,r) = 
^r(m — r) + min(r, m — r). LetFi : D — > W and F 2 : D — > IR m ~ r be the 
functions made up of the first r and the last m — r components of IFF. The 



function F : D — > 
Fix) = 



defined as 




F(x) 
GFAx) 



with G = R\ 2 R^ 



(3.2) 



is "close" to the function F, as shown in the following theorem. 



Theorem 3.1. Let A C Q e f]D be a clos ed convex set containing the origin 
and let F be the function defined by (3.2); for each x G A we have 

\\F(x) - U t F(x) || 2 < q{m, r)5 + 0(5 2 ) 

Moreover, if there exists x* G A such that F(x*) = 0, then 

\\F(x)-Tl t F(x)\\ 2 = 0(5 2 ) 



Proof. W.l.o.g. in (3.1) we may assume II equals to the identity matrix; us- 
ing (3.1) and an obvious partitioning of Q = (Q\ \ Q 2 ) the matrix (Ji?(0)|F(0)) 
can be rewritten as 



(MO) I F(0)) 



^(0) F(0) 
J F2 {0) F 2 {0) 



R n Q[ 

R12Q1 ~l~ R22Q2 



(3.3) 



We denote by d(x) = (F 2 - F 2 ){x) = {GF 1 - F 2 )(x) and so from (|3_2) we 
have that \\F(x)-F(x)\\ 2 = \\F 2 (x) - F 2 (x)\\ 2 = \\d(x)\\ 2 . Further from [3^ 



(GJ Fl (0)-J^(0)\GF 1 (0)-F 2 (0)) 
G(J Fl (0) I F x (0)) - (J Fa (0) |F 2 (0)) 



(J rf (0) I d(0)) = 

= GR\iQ\ — R\ 2 Q\ — R 22 Q 2 = —R 22 Q 
therefore ||J<z(0)||2 < q(m,r)5 and ||d(0))|| 2 < q(m,r)5. 
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We consider the Taylor expansion of d(x) at the origin; since d G C 2 (A) 
and A C Q e n D is a closed and convex set, we have d(x) = d(0) + Jd(0)x + 
0(e 2 ), for each x G A. Since \\x\\^ < y/ne and 5 > e we have 

\\d(x)\\ 2 < \\d(o)\\ 2 + \\j d (o)\\ 2 \\x\\ 2 + o(e 2 ) 

< \\d{0)\\ 2 + q{m,r)^i5e + O{e 2 ) 

< \\d(0)h + q(m, r)^i5 2 + 0(5 2 ) = \\d(0)\\ 2 + 0(6 2 ) 

From ||g?(0))|| 2 < q(m,r)5 the first inequality follows. 

Now, let x* G A such that F(x*) = 0; evaluating the previous Taylor 
expansion of d(x) at x* we get = d(0) + Jd(0)x* + 0(e 2 ), and so d(0) = 
-J d (0)x* + 0(e 2 ). Since \\x*\\ 2 < y/ne, it follows that 

||d(0) || 2 < q(m, r)^i5e + 0{e 2 ) < q{m, r)^i5 2 + 0{5 2 ) = 0{5 2 ) 

and this concludes the proof. □ 



From Theorem 3.1, the partition given in (3.1) can be used in the Root 
Finding Algorithm to derive the new function F. 

Since, by construction, the last m — r equations of the nonlinear system 
F = linearly depend on the first r equations, in order to approximate the 
solution of F = we only have to consider the nonlinear system Fi = 0. 



From formula Q, since (J Fl (0) | Fi(0)) = R^Ql it follows that 

a r (( J>(0) | F(0))) 



a r ((J Fl (0) | Fi(0))) = a^Ru) > 



q(m, r) 



and so (Jp 1 (0) | F x (0)) is well-conditioned if the r-th singular value of the ma- 
trix (Jf(0) I F(0)) is large enough. Nevertheless this fact does not guarantee 
that the matrix Jfi(0) is well conditioned too. For this reason, the Root 
Finding Algorithm, after the construction the function Fx, checks whether 
the numerical rank of Jjr^O) is equal to r. If rank^ (Jf x (0)) < r then Jfi(0) 



is ill-conditioned (see Remark 2.5) and, in order to avoid numerical instabil- 
ities in the computations, no iterations are performed. Otherwise, that is if 
rank^ (J Fl (0)) = r the matrix Jf\(0) is well-conditioned; in this case, start- 
ing from the zero vector, the Root Finding Algorithm processes the system 
F\ = by mean the NFA with an additional check, at each step, on the 
conditioning of the Jacobian matrix, in order to deal with well conditioned 
linear systems. When the loop stops, the algorithm returns the final iterate 
if it satisfies the constraints or, otherwise, a warning message. 



12 



Algorithm 3.2. (The Root Finding Algorithm - RFA) 

Let D C R n be a set containing the origin and F : D — >• M. m with m < n be 
a nonlinear function of class C 2 (D); let w 1 be a fixed threshold, e > 
and Q £ = {x G IR n : ||a;||oo < e}; let 5 > £ of the same order of magnitude 
as £ and k > 1. Consider the following sequence of steps. 

RF1 Compute the numerical (5, /c)-rank of (Jp(0) | F(0))* and a permutation 
matrix II G Mat m (IR) as given by Theorem 2.8 w.r.t. the numerical 
rankr = rank 5 , fc ((J F (0) | ^(0))*). Partition n as (III | n 2 ) with IIi G 
Mat mxr (M); let F x = U\F. 

RF2 Let x = 0, h = (1, . . . , 1)*. 

While (|| ^ || 2 > and (rank,5 )fc ( J Fl (x)) = r) 

• Compute the minimal 2-norm solution h of Jp^x^jh = —Fi(x). 

• Let x = x + h. 

RF3 If x G H -D return x; otherwise return x =NULL. 

Step RF1 requires the determination of the permutation matrix IT; the 



existence of II is guaranteed by Theorem 2.8, its computation can be per- 
formed using a Rank- Revealing QR Factorization [H]. We refer to [3], jUJ for 
a detailed description of algorithms for computing such a decomposition. 

Note that, if Jfi(0) is well-conditioned then there exists a closed neigh- 
bourough A of the origin contained in Q £ n D such that Jf^x) is well- 
conditioned for each x G A . In fact, since F\ G C 2 (D), then Jp^x) is a Lip- 



schitz function (with constant 7) on Ao and so, from Proposition 2.7, for each 
x G A we have that |o>( Jp^x)) — a r (Jp 1 (0))| < || Jf^x) — J^i (0) H2 < Tll^lh 
which implies cr^Jp^x)) > cr^Jp^O))— j\\x\\z- Since a r (JF 1 (0)) > if | |xet|| 2 
is small enough, that is if x belongs to a suitable neighbourough A of the ori- 
gin, then J Fl (x) has full numerical (5, fc)-rank on A and so the RFA executes 
some iterations at step RF2 for approximating the solution of iq = 0. 

We analyze the output of the RFA; if the output x 7^ NULL then the 
peculiarities of x give us different information on the underdetermined non- 
linear system F = 0. If x = is returned, then two possible cases can occur: 
either a single iteration or no iteration of the loop of step RF2 has been 
performed. In the former case h is the zero vector: it means that iq(x) = 
and so F 
Theorem 



x) = 0, that is the exact solution of F = has been found. From 
it follows that ||F(x)|| 2 = \\F(x) - F(x)\\ 2 < q(m,r)5. In the 



3.1 



latter case no iteration has been executed since rank^J^O)) < r, though 
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rank,5 jfc ( Jp(0) | F(0)) = r; in this situation h — (1 ... 1)* and, if some simple 



hypotheses on a closed convex set A C Q e n D are satisfied, Theorem |3.3 
shows that there are no solutions of F = in A. 

Theorem 3.3. Let A C Q £ DD be a closed convex set containing the origin 
and let 7 be the Lipschitz constant of J Fl on A. Suppose that the matrix 
(J Fl (0) I Fi(Q)) has full numerical (5, k)-rank r with k > 1 + (s/n + wj)S. If 
rank^ (J Fl (0)) < r then there are no solutions of F = in A. 

Proof. Note that, since F G C 2 (A) and A is a closed convex set, J F is a 
Lipschitz function on A. Let r\ = rank^ («7fi(0)) and A G Mat rXn (R) be 
such that rank(A) = r± and 1 1,7^(0) — A\\ 2 < 5: the existence of such a 



matrix A follows from Lemma |2.6[ Since A has exact rank ri there exist a 
permutation matrix LT G Mat r (R) and a matrix W G Mat( r _ ri ) xri (M) such 
that 

Ax \ }n 



HA = \ rxr a i an d rank(Ai) = n 

We can express the matrix HJ Fl (0) as follows 

nj«(0) = ( ^+ + ^ 21 ) with B = ( |; ) and < S 

Suppose for a contradiction that there exists x* G A such that F(x*) = 0, 
and so Fi(x*) = 0. From Taylor's theorem applied to Fi(x) we know that 
there exists a point z on the line connecting to a;*, such that = F\(x*) = 
Fi(0) + J Fl (z)x*, and so 

= nFi(o) + tu Fi (z)x* = nFi(o) + nj Fl (o)x* + n(j Fl (z) - J Fl {o))x* 

Using an obvious partitioning (Fn,Fi 2 Y and (Ei 2 ,E 22 Y of the matrices 
nFi(O) and n(J^(z) - J Fl (0)) we obtain 



and therefore 



A lX * = -F u - (E n + E 12 )x* 

F12 = W(Fn + (E n + E 12 )x*)-(E 21 + E 22 )x* 
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It follows that II (Jp 1 (0) | -Fi(O)) is equal to the matrix 

A\ + E\\ F\\ 
WA 1 + E 21 W (F n + (E u + E 12 )x*) - (E 21 + E 22 )x* 



Let B G Mat rx „(R) such that 
UB = 



A 1 F u + (E u + E l2 )x* 

WA X W (F n + (E u + E 12 )x*] 

since || • H2 is invariant under permutation and x*, z G A we have 

\\B - (J Fl (0) \ F 1 (0))\\ 2 = \\(-E\Ex* + Tl(J Fl (z)-J Fl (0))x*)\\ 2 

< \\E\\ 2 (1 + ||x*|| 2 ) + \\J Fl (z) - J Fl (0)h\\x*h 

< 5(1 + y/ne) +^ne 2 < 5 + (y/n + wy)5 2 

From rank(A!) = r\ we have that rank(5) = r\, and so <J r (B) = 0. From 
Proposition |Z7| it follows that o, r (J Fl (0) \ Fi(0)) < 5 + (y/n + n^)5 2 . Since 
from the hypothesis a r (J Fl (0) |Fi(0)) > k5 and k > 1 + (y/n + n^)8 we obtain 
a contradiction. □ 

If the RFA returns a non zero vector x G Q £ fl D then F(x) assumes 
small values; in particular, ||F(a;)||2 is bounded by a function depending on 
the data tolerance as the following theorem shows. 

Theorem 3.4. Let x be the output of the RFA applied to F G C 2 (D) and H 
be the closed convex hull of the sequence of points computed by the RFA. If 
H C Q £ H D then 

\\F(x)\\ 2 = 0(5 2 ) and \\F(x) || 2 < q(m, r)6 + 0(5 2 ) 

Further, let A be a closed convex set such that H C A C Q E n D; if there 
exists x* G A such that F(x*) = then ||F(x)|| 2 = 0(5 2 ). 

Proof. Let h be the displacement computed in the last iteration of the RFA, 
that is let h satisfies J Fl (x — h)h = —Fi(x — h) and consequently Jp(x — 
h)h = —F(x — h). Since x — h G H, F G C 2 (H) and if is a closed convex 
set, we express F(x) using the Taylor expansion of F centered at x — h 
we have F(x) = F(x - h) + J p (x - h)h + 0(\\h\\ 2 ) = 0(\\h\\ 2 ) and, from 
\\h\\ 2 < v^e < y/n~5, it follows that ||F(z)|| 2 = 0(5 2 ). 

Since \\F(x)\\ 2 = ||nF(x)|| 2 < \\F(x) - ILF(x)\\ 2 + \\F(x)\\ 2 , from Theo- 
rem 3.1 we have ||F(x)|| 2 < q(m, r)S + 0(S 2 ) and, if there exists x* G A such 
that F(x*) = then \\F(x)\\ 2 = 0(5 2 ). □ 
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4 The Low-degree Polynomial Algorithm 



Given a set X e of s distinct empirical points of M. n in this paper we address 
the problem to determine a polynomial of P = WL[xi, . . . , x n ] whose degree is 
bounded by the minimal degree of the elements of the vanishing ideal X(X) 
and whose affine variety exactly contains an admissible perturbation X* of X E . 

In this section we present an algorithm that computes an admissible per- 
turbation X of X £ and a polynomial /, whose degree is bounded by the 
minimal degree of the elements of X(X); since / has low degree and, by con- 
struction, assumes small values at X, it is a good candidate to be the solution 
to the addressed problem. In fact, if some simple additional hypotheses are 
satisfied, / turns out to be the polynomial we are looking for, as in Theo- 
rem 4J5 we prove that there exists an admissible perturbation X*, slightly 
differing from X, contained in the affine variety 2(f). 

As already pointed out in the Introduction a straightforward application 
of the Buchberger-Moller (BM) algorithm to empirical data provides unre- 
liable results, because the problem of deciding whether a polynomial van- 
ishes at the input points is very sensitive to perturbations. Nevertheless the 
strategy of the BM algorithm may be adapted to limited precision data; for 
instance the NBM and SOI algorithms (j2], [ID]) employ such an approach to 
compute polynomials which almost vanish at the input points. Also in this 
context we present a modified version of the BM algorithm. 

The core of the BM algorithm is the stepwise construction of the mono- 
mial basis O of the quotient ring P/X(X) viewed as an IR-vector space. 
Initially O comprises just the power product 1; at the generic iteration 
O = {ti, . . . , ty) and, fixing a term ordering r, the smallest t > T t i: i = 1 . . . v, 
is considered. If the evaluation vector t(X) is linearly dependent on the 
columns of the matrix Mq (X) , that is if 

p(X) = (4.1) 

where p(X) = t(X) — Mo(X)a(X) is the residual of the least squares problem 
Mo(X)ot(X) = t(X), then the polynomial t — J2i=i a iOQ^i * s formed and 
added to a basis, otherwise t is added to O. 

We adapt the generic iteration of the BM algorithm to our case by re- 



placing check (4.1) by 



p(X*) = for some admissible perturbation X* of X £ (4-2) 
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If problem (4.2) admits solution then the polynomial f* = t — Yli=i 

is formed; /* and X* are a solution we were looking for (by construction /* 

vanishes at X*). Otherwise t is added to O. 



Analogously to the strategy used in [2], in order to make problem (4.2) 
effectively solvable we express any admissible perturbation of X e as a function 
of sn error variables e = (en, . . . , e sl , e 12 , . . . , e s2 , • • • , e lnj . . . , e sn ), where 
each variable e k j represents the perturbation in the j-th coordinate of the 
specified value p k G X. Specifically, a generic admissible perturbation X 
of X £ can be expressed as X = X(e) = {pi(ei), . . . ,p s (e s )} where for each 
k = l...s 

Vki^k) = (Pki + e kl , p k2 + e fc2 , • • -Pkn + e kn ) and He^oo < e 

Using this formalization, X = X(e) is a function in the variables e with 
admissible domain Q £ = {e G lR sn : HeU^ < e}. It follows that ct(X), p(X), 
and Mq(H) are functions of e, simply denoted by a(e), p(e) and Mo(e); in 
particular p : W n — > W is defined by 

p(e)=t(e)-M (e)a(e) (4.3) 



where a(e) = Mo(e)H(e). It follows that problem (4.2) is equivalent to 
determine whether 

p(e) = for some e e Q £ nD (4.4) 
where Do is the domain of p, that is Do = {e G lR sn : Mp(e) is full rank}. 



In [10] it is proved that the problem (4.4) has no solution if 



|p(0)| > e\I - M o (0)M o (0y\ ^ \d k t(0) - M SfcO (0)a(0)| + 0(s 2 ) (4.5) 



k=l 



where the absolute value of a matrix is the matrix consisting of the absolute 
values of its elements and the lower bound means that the relation holds com- 



ponentwise. In our algorithm we employ such a result to study problem (4.4 ): 



we solve the underdetermined nonlinear system (4.4) only if condition (4.5) 
is not satisfied. Unfortunately, finding an exact solution of problem (4.4) can 
be computationally very hard: in fact the analytic expression of p(e) can be 
very hard to obtain and even when an explicit expression of p(e) is available, 
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an exact solution e* of p(e) = can be computed only in very few cases 
(specifically when n and s are small), while in general only an approxima- 
tion of e* can be effectively determined. In our algorithm we approximate 



the solution of problem (4.4) using the results of Section 3 in particular, we 



design the The Low-degree Polynomial Algorithm following the frame of the 



BM algorithm and solving the nonlinear systems (4.4) using the RFA 



Algorithm 4.1. (The Low-degree Polynomial Algorithm - LPA) 

Let X e = {pi, . . . ,p £ s } be a finite set of distinct empirical points of M. n , k > 1, 
5 > e of the same magnitude as e and r be a term ordering. 
Consider the following sequence of steps. 

LP1 Start with t = 1 and the lists O = [1] and L = [ ]. 

LP2 Add to L those elements of {xit, . . . ,x n t} which are not multiples of 
elements of L. Let t = min r (L) and delete it from L. 



LP3 Choose e = and consider the residual function p defined in (4.3) 



LP4 If p(0) satisfies condition (4.5) then 

• add t to the list O and continue with step LP2; 
else 

• apply the RFA to p on Q £ fl Do with thresholds k and S. 

LP5 If the RFA returns NULL or e = and \\h\\ ^ then 

• add t to the list O and continue with step LP2; 
else 

• compute / = t — Yliti&o return e, /, O and stop. 

We make a few observations on the LPA. At each iteration the domain Do 
of the residual function p contains the origin. In fact since the residuals 
computed at the previous iterations do not vanish at zero the matrix Mo{0) 
is full rank. Further, at step LP4 each iteration of the RFA requires the 
evaluation of p and J p at a given point e. The value of p(e) is simply 
obtained by solving the numerical least squares problem Mo(e)a(e) = t(e) 



and computing its residual. Proposition 6.1 illustrates how to compute each 
column of the Jacobian matrix J p (e) without passing through the explicit 
expression of the function p. 

Theorem 4.2. The Low-degree Polynomial Algorithm stops after finitely 
many steps and returns e, f , O such that eeQ £ n Do- 
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Let p be the residual function computed during the last iteration of the 
LPA and H be the closed convex hull of the sequence of points computed by 
the RFA applied to p. If H C Q £ n Do then 

\\f(X(e))h<q(s,r)S + 0(5 2 ) 

where r = rank^ J p (0) | p(0)) and q(s, r) = a/ r(s — r) + min(r, s — r). 

Further, let A be a closed convex set such that H C A C Q £ fi Do; if 
there exists e* G A such that p(e*) = then ||/(X(e))|| 2 = 0(5 2 ) 

Proof. At each iteration either the algorithm computes the polynomial / 
and thus it stops, or a term t is added to O. We observe that this latter 
instruction can be executed at most s — 1 times; in fact when O contains s 
terms, that is when Mq(0) becomes a square matrix, the residual vector p(0) 
is zero, and consequently e = 0, p(e) = 0, and a polynomial / is computed. 
From the check of step LP 5 obviously follows that e G Q e PI Do- 

The residual function p e C°°(Do) since each of its components is a 
rational function defined in Do- The thesis follows from Theorem |3 . 4| applied 
to p and the equality /(X(e)) = p(e). □ 

In the following theorem, under additional hypotheses on the polynomial 
/, output of the LPA, we prove that there exists an admissible perturba- 
tion X(e*) of X £ contained in Z(f); further, we provide an estimation of the 
distance between X(e) and X(e*) which is usually much smaller than e. 

Theorem 4.3. Let (e,f,0) be the output of the LPA applied to X e . For 
each % = 1 . . . s let fi : M. n — > R be the function defined by /i(ej) = f{pi{<£i)); 
let Ti and Ri be real positive numbers such that 

f 11^(^)11 

Ti = e — 1 1 e ? ; 1 1 oo and Ri < min < r^, — - — - — 

Ti 



where ji > is a Lipschitz constant of Jf in B{e i) r i ). Let /ij be an upper 

'I 



bound of \[j\.\[2 i n B{e. u Ri). If \\JfX e i)\\2 > and 



l/t(e»)l< (0 ^\ ^ = Xi Vi = l...s (4.6) 

Hi(2 + "fiRiHi) 

then there exists e* such that the set X(e*) is an admissible perturbation 
ofX £ , /(X(e*)) = ; and \\e* — e||oo < maxiZj. 
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Proof. Let dj = <x,(e) and f = t — Y^fyeo^j- We observe that /i(ej) = 
t{Pi{ e i)) + Yltjeo ®jtj(Pi( e i)) * s a polynomial function of and its Jaco- 
bian JfX e i) ^ Mati xn (R) is a Lipschitz function in 5(e«,rj). We prove that 
l|J/i( e «)ll2 > in B(ei,Ri). For each G fi(e i ,i? i ) we have 

I 11^(^)11 2 - \\Jfi(ei)h | < \\Jfi(ei) - J fl (^)h < liW^i - < ^Ri 

and so ||«// 4 (ej)|| 2 > || J^(ej)|| 2 — jiRi > 0. It follows that Jfi( e i) nas constant 
rank 1 in B(ei,Ri); further, (e^ ) 1 1 2 is upper bounded in B{e i: Ri) since 

+ II^UeJL 1 1 



JfMWl \\ J fi( e i)h minB^A) H J /,( e i)ll2 

and min 5(giA) || Jf^ijh > 0. 

Now we apply Kantorovich theorem to /, on B(ei,Ri) with x = e 



2, 



p = 1, r] = Ri and so = {ei}; since Xi = Mi ( 2 +7-fl iMi ) satis fi es "^t^ 1 < 1 
and ,„^ iX i ■. < Ri, and (14.6 ) holds we conclude there exists e* = (e% . . . e* ) 

in B(ei, Ri) such that /i(ef) = 0. Let e* = (e* n . . . e* sl , . . . , e* n . . . e* n ); since 
/(X(e*)) = (/i(e*), . . . , / s (e*))* we have /(X(e*)) = 0. Furthermore, since 
the i-th element of X is pi = Pi(0) and 

\\Pi(e*) -p^OJIloo = lle-Hoo < ||e* - §j|| a + IJejUoo < Ri + INU < e 

we have that X(e*) is an admissible perturbation of X e . Finally, we have 
that ||e* — e||oo = maxj ||e* — ej||oo < max; ||e* — || 2 < maxj/2j < e which 
concludes the proof. □ 



5 Numerical examples 

In this section we present some numerical examples to show the effectiveness 
of the LPA. We implemented the LPA and the RFA using the C++ language, 
the CoCoALib [5], and some routines of GSL - GNU Scientific Library [12]; all 
computations have been performed on an Intel Core 2 Duo processor (at 1.86 
GHz). In all the examples the degree lexicographic term ordering with y < x 
is understood; furthermose the coordinates of the points and the coefficients 
of the polynomials are sometimes displayed as truncated decimals. 

In Example 5.1| the LPA is applied to the set of points of Example 1.1 



created by perturbing by less than 0.1 the coordinates of 10 points located 
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on the parabola g = where g = y 2 — x — 2y + 2. In contrast to the exact 
approach that provides a minimal degree polynomial of degree 4, the LPA 
detects the geometrical configuration given by a parabola. 

Example 5.1. Let X C 1R 2 be the set of points given in Example let 
e = 0.1, 5 = 2e and k = 2. The LPA applied to X £ performs three iterations. 
At the first two iterations the terms t — y and t = x are analyzed: since in 



both cases condition (4.5) is satisfied the set O = {l,y,x} is constructed. 
During the third iteration the term t = y 2 is considered: at step LP4 the 
RFA applied to the current residual function p returns e e Q £ nDc, therefore 
at step LP5 a polynomial / is computed (see Figure [T]), and the algorithm 
stops with the following output: 

e = 10- 1 (0.474, -0.211, -0.146, -0.100, 0.075, -0.105, -0.091, -0.017, 
0.066, 0.054, 0.004, -0.215, 0.285, -0.202, 0.211, 0.378, -0.332, 
-0.065,-0.369,0.305) 

/ = y 2 _ o.9751012065x - 2.0049270587y + 1.9775224038 
O = {l,y,x} 



Figure 2: The varieties Z(f) and 
Figure 1: Variety Z(f) Z{K): zoom in around the point 

(0.95, 1) 

The set X(e) is an admissible perturbation of X £ ; we can consider the 
polynomial / almost vanishing at X(e) since ||/(X(e))|| 2 /||/|| < 2 • 10~ 5 . We 



apply Theorem 4.3 to this example and report the values of Ri, Xi and |/i(ej 



in the following table. Since for each i the inequality (4.6) is satisfied, we 
conclude that there exists e* such that X(e*) is an admissible perturbation 
of X e , /(X(e*)) = and ||e* - < 0.0935. 
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i=l 


i=2 


i=3 


i=4 


i=5 


Ri 


0.0526 


0.0785 


0.0715 


0.0798 


0.0789 


Xi 


0.0243 


0.0354 


0.0325 


0.0360 


0.0356 


\fi(*i)\ 


3.3- 10~ 7 


3.87- 10~ 6 


3.46 • 10~ 6 


4.77- 10~ 6 


1.46 • 10~ 6 




i=6 


i=7 


i=8 


i=9 


i=10 


Ri 


0.0622 


0.0668 


0.0935 


0.0631 


0.0695 


Xi 


0.0285 


0.0305 


0.0416 


0.0289 


0.0316 


\fi(ei)\ 


9.23 • 10" 6 


4.45 • 10" 6 


1.015 • 10" 5 


2.783 • 10" 5 


4.95 • 10" 6 



The NBM algorithm as well as the SOI algorithm (available in [nj) ap- 
plied to X £ detect the same support of /, since they compute the polynomial 
h = y 2 - 1.0041s - 2.0089y + 2.1287 which almost vanishes at X. Neverthe- 
less h is not a solution to our problem: in fact, h = (the dashed line in 
Figure [2]) contains no admissible perturbation of X £ since it does not cross 
the e-neighborhood of the first point of X. 

In the following example a fairly large set of points of M. 5 is considered. 

Example 5.2. Let X C M. 5 be a set of points created by perturbing by less 
than 0.1 the coordinates of 25 points lying on the affine variety Z(g) where 

fi ™ 2 1 7 ,y, ,y, o ,y, 2 10 ,y, i 21 y, 74 ™ I 93 ,y, 36 ™ i 39 

i/ — • i '4 4x Jy 4- x -5 ^-^5 4l x l ' 41 2 4i Jy 3 ~r 41 -^4 4i x 5 T" 41 

Let e = 0.1, 8 = 2e and k = 2; the LPA applied to X £ performs 11 iterations 
in a computational time of 24 seconds. During the 11th iteration the term 
t = x\ is considered and the output e, /, O is returned, where He]^ = 0.0789 

/ = x\ - 0.421x 4 x 5 - 1.994x| - 0.244xi + 0.405x 2 - 1.779x 3 + 2.328x 4 
-0.685x5 + 1-124 + 10- 2 • (0.056xix 5 - 2.549x 2 x 5 + 1.679x 3 x 5 ) 

O = {l,X5,X4,X3,X2,Xi,x|,X4X5,X 3 X5,X2X5,XiX 5 } 

The set X(e) is an admissible perturbation of X e ; we can consider the poly- 
nomial / almost vanishing at X(e) since ||/(X(e))|| 2 /||/|| ~ 3.42 ■ 10" 4 . We 



apply Theorem 4.3 to this example; since for each i the inequality (4.6) is 
satisfied, we conclude that there exists e* such that X(e*) is an admissible 
perturbation of X e , /(X(e*)) = and ||e* - < 0.097. 

We observe that the minimum of the degrees of the polynomials of X(X) is 
3, and this does not suggest that the points of X lie close to the variety 2(g). 
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On the contrary, since the coefficients of / differ only slightly from the cor- 
responding coefficients of g, the LPA allows us to recover the "approximate" 
geometrical configuration of X. 

The following example (taken from [7]) is different from the standard 
cases deriving from the analysis of real-world phenomena. It points out that 
the LPA has various possibilities of applications; in this case it is employed 
to obtain the numerical implicitization of a parametric curve. 

Example 5.3. We consider the parametric equations for a Bezier curve: 

x{t) = At(2t 5 - 3t 4 + 8t 2 + 6t + 3)/(t 6 - 3t 5 + 3t 4 + 3t 2 + 3t + 1) . . 
y{t) = 6t(4t 4 + 9t 3 - 9t 2 -9t + 5)/(t 6 - 3t 5 + 3t 4 + 3t 2 + 3t + 1) [ } 

whose implicit equation is g = where 

a = 2 T 2„ 28 T7/ 2 i 224 - 3 15712 2 56 i 848 „ 2 

y .x, 12 69 J ' w 423 y ~ 34263 y 1269 1269 y ' 3807 y 

I 44480 _ 17792 

' 1269 1269 * 

We consider the set of points 

X = {(0,0), (-1.3581, -4.7661), (2.0956, 2.0315), (4.6884, -0.3349), 
(-2.7205, -11.6848), (-7.2835, -40.9773), (6.7793, -1.4114), 
(8.6024, 1.4575), (10.5250, 8.8937), (12.6213, 19.7217)} 



created by evaluating the parametric equations (5.1) at 10 random values of 
the parameter t G (—1, 2), and rounding off up to 10~ 4 . 

Let e = 10~ 4 , 5 = 2e and k = 2; the LPA applied to X e performs 9 it- 



erations. During the first 8 iterations condition (4.5) is satisfied, and the 
set O = {l,y,x,y 2 ,xy,x 2 ,y 3 ,xy 2 ,x 2 y} is constructed. During the 9-th iter- 
ation the term t = x 3 is considered: at step LP4 the RFA applied to the 
current residual function p returns §i G Q £ fl Dq, therefore at step LP5 a 
polynomial fi is computed, and the algorithm stops with the output: 

e x = 10- 6 - (8.13,-1.99,0.92,-4.90, 1.40,-7.63,-3.30,23.79,-24.90, 
8.50, -3.25, 0.36, -0.13, 1.11, 4.73, -7.19, 9.02, -8.41, 5.32, -1.54) 

f x = x 3 - 0.00159x 2 ?/ - O.O66I82?/ 2 + 0.00654?/ 3 - 12.3814a; 2 - 0M396xy 
+0.22269?/ 2 + 35.0514a; - 14.0206?/ - 0.00033 

O = {l,y,x,y 2 ,xy,x 2 ,y 3 ,xy 2 ,x 2 y} 
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Note that the coefficients of fi and g are very close to each other; nev- 
ertheless, an even better result can be obtained by taking into account the 
exact information that the Bezier curve g = vanishes at the origin. To 
this aim we impose that the constant term of the polynomial computed by 
the LPA is zero; under this assumption and using the above values for the 
parameters e, 5, k, the LPA applied to X £ stops with e 2 , f 2 , O, where 

e 2 = 1(T 6 ■ (0, -2.09, 0.96, -5.13, 1.46, -8.00, -3.45, 24.92, -26.09, 8.90, 
0, 0.37, -0.14, 1.16, 4.95, -7.53, 9.44, -8.81, 5.57, -1.61) 

f 2 = x 3 - 0.00158x 2 ?/ - 0.06619x?/ 2 + 0.00653?/ 3 - 12.3814x 2 - 0.04407xy 
+0.22271?/ 2 + 35.0513x - 14.0205?/ 

We observe that f 2 is better than fi because its coefficients provide an even 
closer approximation of the coefficients of g, and its support and the support 
of g are the same. 



6 Appendix 

The following proposition illustrates how to compute each column of the 
Jacobian matrix J p (e), that is the vector dp(e)/deki, without passing through 
the explicit expression of the function p. 



Proposition 6.1. Let O — {t\ . . . t m } be a finite set of terms, let t be a term, 



t ^ O, and let p be given by (4-3); then 
dp(e 



de k ,i 



I — Mo(e)Mo(e) 



dt(e) dM (e) 



de 



hi 



de 



ki 



deu 



-P(e) 



where I G Mat s (JR) is the identity matrix, M J (e) = (M|,(e)M (e)) 1 M^(e), 
dt(e)/deki is the zero vector except the k-th coordinate which is equal to 
dit{pk{ G k)) and dMo(e)/deki is the zero matrix except the k-th row which is 
equal to (9^i(p fc (e fe )) . . . dit m (p k (e k ))) . 

Proof. From the equality p(e) = t(e) — Mo(e)a(e), deriving matrices and 
vectors componentwise w.r.t. e^, we obtain 



dp(e) _ dt(e) 
de k i de 



dM (e) da(e) 
-a(e) — %le 



ki 



de 



ki 



ki 



(6.1) 
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Deriving componentwise each side of Mg(e)M c (e)a(e) = Mp(e)f(e) w.r.t. e ki 
we have 



da(e) 



(MS(e)M o( e)) -^0) + M> o( e) ^ _ _^Li a(e ) j 



and the thesis follows by substituting the expressions of da(e)/de k i in ( |6.1| ). 
|M and 



For describing and ° ^ we consider the derivative of a generic 



term q = nf =1 xf\ Since g(e) = (n? =1 (pii + e u f% • • • , n^ifr* + e si ) ft ) the 
vector is the zero vector except its fc-th coordinate which is equal to 
diq(p k (e k )). It follows that only the fc-th coordinate of the vector and 
the fc-th row of the matrix 9K q°^ can be different to zero: they are equal to 
dit{p k {e k )) and to (9<*i(pfc(e fc )) . . . dit m {p k {e k ))) respectively. □ 
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